Assessment of density functional approximations for the hemibonded structure of 

water dimer radical cation 



(N 

o 

(N 

< 



43 

Oh. 



Piin-Ruey Pan, 1 You-Sheng Lin, 1 Ming-Kang Tsai, 2 Jer-Lai Kuo, 3 'B and Jeng-Da Chai 1,4 '0 

1 Department of Physics, National Taiwan University, Taipei 10617, Taiwan 

2 Department of Chemistry, National Taiwan Normal University, Taipei 11677, Taiwan 

3 Institute of Atomic and Molecular Sciences, Academia Sinica, Taipei 10617, Taiwan 

4 Center for Theoretical Sciences and Center for Quantum Science and Engineering, 

National Taiwan University, Taipei 10617, Taiwan 
(Dated: March 2, 2013) 

Due to the severe self-interaction errors associated with some density functional approximations, 
conventional density functionals often fail to dissociate the hemibonded structure of water dimer 
radical cation (H20)2 + into the correct fragments: H2O and EbO" 1 . Consequently, the binding 
energy of the hemibonded structure (H20)2 + is not well-defined. For a comprehensive comparison 
of different functionals for this system, we propose three criteria: (i) The binding energies, (ii) the 
relative energies between the conformers of the water dimer radical cation, and (iii) the dissociation 
curves predicted by different functionals. The long-range corrected (LC) double-hybrid functional, 
cjB97X-2(LP) [J.-D. Chai and M. Head-Gordon, J. Chem. Phys., 2009, 131, 174105.], is shown to 
perform reasonably well based on these three criteria. Reasons that LC hybrid functionals generally 
work better than conventional density functionals for hemibonded systems are also explained in this 
work. 
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I. INTRODUCTION 

Water can be decomposed when it is exposed to high- 
energy flux. The products of water radiolysis may con- 
tain various radical species, e.g. hydrogen atoms (H), 
hydroxide radicals (OH), oxygen anions (O - ), and wa- 
ter cations (H2 + ), depending on the radiation infras- 
tructure setup. For example the overall decomposition 
scheme activated by /3 particles has been outlined by 
Garrett et at in 2005 [l| where three main channels of 
decomposition were listed. The cationic channel leads 
to the formation of ionized water living in about sev- 
eral tens femtoseconds and hydrated electron, followed 
by the generation of hydronium (HaO + ) and OH radicals 
through proton transfer process @, The energized- 
neutral and anionic channels could result in the cleavage 
the oxygen- hydrogen chemical bonds to produce hydro- 
gen and oxygen derivatives, i.e. H, H~, H 2 , O, O - , OH~ 
etc. Subsequent chemical reactions can progress further 
up to the desorption of stable gas molecules H2 and O2 
being driven by those reactive radical species The 
cationic channel is therefore particularly interesting due 
to its dominant products — OH radicals and solvated 
electrons. 

The smallest system to understand the chemical dy- 
namics of ionized water is the water dimer radical cation 
(H20)2 + , and it has been approached by several exper- 
imental studies in the past. Angel and Stace reported 
the predominant HaO + -OH central core by the collision- 
induced fragmentation experiment [|| against the ear- 
lier theoretical assignment of charge-resonance hydrazine 
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structure [5j. Dong et al. observed a weak signal corre- 
sponding the formation of (H 2 0) 2 + near the low-mass 
side of (H 2 0)2H + using 26.5 eV soft X-ray laser [6|. Gar- 
denier, Johnson, and McCoy reported the argon-tagged 
predissociation infrared spectra of (H2 0)2 + and assigned 
its structural pattern as a charge-localized H30 + -OH 
complex Q. Recently, Fujii's group reported the infrared 
spectroscopic observations of larger (H2 0)„ + clusters, 
n = 3 — 11 8], where the OH radical vibrational signal 
was clearly identified for n ^ 5 clusters, but vibrational 
signature of OH radical become inseparable due to the 
overlap with H-bonded OH stretch in n > 6 . As be- 
ing evidenced in the earlier studies 0, HI , theoretical in- 
vestigations such as ab initio electronic structure theory 
and density functional theory (DFT) play an important 
role in understanding the infrared spectroscopic features 
of the ionized water clusters. Because high-level ab ini- 
tio calculations are computationally prohibited for larger 
ionized water clusters, e.g. fully solvated cationic moiety, 
a reliable DFT method is necessary. 

In earlier theoretical reports, two minimum structures 
of the water dimer radical cation have been identified: the 
proton transferred structure and the hemibonded struc- 
ture, as shown in Fig.[Tl[9rll2l|. The previous DFT calcu- 
lations have shown that many exchange-correlation (XC) 
functionals fail to predict reasonable results 0-[ll| given 
rise to the presence of hemibonding interaction. The 
hemibonding interaction, could be theoretically located 
in (H2 0)„ + systems, is notorious for the serious self- 
interaction errors (SIEs) associated with some density 
functional approximations. Both local density approxi- 
mation (LDA) and generalized gradient approximations 
(GGAs) were reported to contain non-negligible amount 
of SIEs for describing the hemibonded structure [9|, llOj . 
It has been suggested to adopt hybrid functionals with 
larger fractions of the exact Hartree-Fock (HF) exchange 



2 



for more accurate results in the hemibonded structure 
9, 10]. However, as the SIEs of functionals become larger 
at the dissociation limit, these suggested functionals can 
yield spurious barriers on their dissociation curves [To| . 
which can lead to unphysical results in molecular dynam- 
ics simulations. 

Clearly, more stringent criteria for choosing suitable 
functionals are needed. In this work, we propose three 
different criteria for a comprehensive comparison of dif- 
ferent functionals for this system. 



(«) 
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(b) (c) 



FIG. 1. (a) The proton transferred structure, (b) the hemi- 
bonded structure, and (c) the transition state between the 
structures of (a) and (b). 



All of the calculations are performed with the devel- 
opment version of Q-Chem 3.2 [32]. As the basis set 
superposition error (BSSE) for the ionized water dimcr 
has been shown to be insignificant (if the diffuse basis 
functions are adopted) [9j, 111], we do not perform BSSE 
correction throughout this paper. 



III. RESULTS AND DISCUSSION 

The ZPE corrected binding energies and relative ener- 
gies of the water dimer cation calculated by various XC 
functionals are shown in Table Q] and Table [Til respec- 
tively. The calculated dissociation curves for the hemi- 
bonded structure are shown in Fig. [2] A summary of the 
results based on these three different criteria is shown in 
Table HVl The notation used for characterizing statisti- 
cal errors is as follows: Mean signed errors (MSEs), mean 
absolute errors (MAEs) and root-mean-square (RMS) er- 
rors. 



A. Criterion I: Binding Energies 



II. COMPUTATIONAL METHODS 

Calculations are performed on the optimized geome- 
tries of the two structures of the water dimer cation 
and the transition state between them, optimized with 
the ab initio MP2 theor y I13l| and various XC function- 
als involving BLYP [U, [T^TpBE 0, and M06L [nj], 
which are pure density functionals (i.e. the fraction of 
HF exchange a nF = 0.00), B97 [3 with ct H F = 0.19, 

20 with 



22] with 



B3LYP [3 [H [Tlwith a HF = 0.20, PBE0 
«hf = 0.25, M06[2l| with a H F = 0.27, M05 
a H F = 0.28, BH&HLYP || with a HF = 0.50, M06-2X 
[Hi with a H F = 0.54, M05-2X Q with a H F = 0.56, 
M06HF [H] with a H F = 1-00, the wB97 series (wB97 
[26|, wB97X [26|, wB97X-D g^, and wB97X-2(LP) [H), 
which are long-range corrected (LC) hybrid functionals 
(i.e. the fraction of HF exchange depends on the interelec- 
tronic distance [29l|). and the double- hybrid functional 
B2PLYP [H with a H F = 0.53. The DFT and MP2 
calculations are performed with the 6-311+- t-G(3df,3pd) 
basis set, where the reference values of binding energy are 
obtained from Ref. [l2j. For efficiency, the resolution-of- 
identity (RI) approximation [3l| is used for calculations 
with the MP2 correlation (using sufficiently large auxil- 
iary basis sets). 

The CCSD(T) dissociation curves are calculated on 
the fixed monomer geometries (using the CCSD(T) op- 
timized geometry of Ref. [ll|), with the aug-cc-pVTZ 
basis set. Note that although the ZPE corrected energy 
of the proton transferred structure (or, referred to as the 
Ion structure in Ref. [Hj].) is inconsistent with Ref. [13 . 
However, adopting the geometries obtained from Ref. [ll| 
yields results that are consistent with Ref. 12 1. 



TableUshows the binding energies of the two structures 
of water dimer radical cation and the transition state be- 
tween them. The reference data obtained from Ref. [l2[ 
are based on the CCSD(T) calculations with the aug-cc- 
pVQZ basis set. Since the errors of XC functionals for the 
hemibonded structure are much larger than those for the 
proton transferred structure, we focus our discussion on 
the hemibonded structure. From Table U the function- 
als with MAE less than 2.5 kcal/mol are wB97X-2(LP), 
M06HF, and BH&HLYP. Table [Q also confirms the trend 
that has been mentioned previously: Functionals with 
larger fractions of HF exchange give more accurate re- 
sults for the hemibonded structure. To give reasonable 
results for the hemibonded structure, the Ohf of a global 
hybrid functional should be at least larger than 0.43, as 
observed in Ref. [ill for tne MPW1K functional. Al- 
though M06HF, containing a full HF exchange, gives 
a small error of the hemibonded structure, it yields a 
large error for the proton transferred structure (due to 
the incomplete cancelation of errors between the exact 
exchange and semilocal correlation), as shown in Table 
|TJ Although functionals with qhf larger than 0.43 have 
been suggested, functionals with ohf > 0.5 are not al- 
ways reliable, which can be observed from the errors of 
the hemibonded structure calculated by B2PLYP, M05- 
2X and M06-2X. But this trend still holds: the results of 
the M05-2X and M06-2X are much better than those of 
the M05 and M06. This means that although the energy 
of the hemibonded structure is sensitive to the «hf values 
in XC functionals, it may also be affected by the asso- 
ciated density functional approximations (DFAs). Also 
note that some GGA functionals, such as BLYP and 
PBE, cannot predict that the transition state between 
the two structures of the water dimer radical cation. 
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TABLE I. Binding energies (in kcal/mol) of the ionized water dimer. 



Proton transferred structure Transition state 



Hemibonded structure 



A jT 4- 1-. ^ J 

Method 


OfHF 


hj 


Error 


hj 


Error 


ri 


Error 


MSE MAE RMS 




0.00 


-45.62 


2.10 






-52.89 


18.17 








PBE 


0.00 


-47.37 


3.85 






-53.93 


19.21 








M06L 


0.00 


-45.00 


1.48 


-40.85 


12.45 


-48.24 


13.52 


Q 1 Pi 

y. 10 


y. 10 


1U.00 


B97 


0.19 


-45.33 


1.81 


-38.39 


9.99 


-45.88 


11.16 


I .00 


7 

I .00 


ft 71 

O. ( 1 


B3LYP 


0.20 


-45.94 


2.42 


-38.43 


10.02 


-45.67 


10.95 


i .oU 


i .oU 


o.Oo 


PBEO 


0.25 


-46.61 


3.09 


-36.75 


8.35 


-43.95 


9.23 


O.oy 


o.oy 


7 A(\ 


M06 


0.27 


-45.83 


2.31 


-36.13 


7.73 


-42.43 


7.71 


o.yz 


o.yz 


(\ AA 


M05 


0.28 


-45.35 


1.83 


-35.46 


7.06 


-41.31 


6.59 


0. 10 


0. 10 


O.Oo 


TJTTP TTT \/"T> 


0.50 


-45.97 


2.45 


-29.55 


1.15 


-35.11 


0.39 


1.00 


l.OO 


1 

l.Oo 


B2PLYP 


0.53 


-45.04 


1.52 


-32.58 


4.18 


-40.96 


6.24 


O.yo 


o.yo 


A AO 


M06-2X 


0.54 


-47.05 


3.53 


-32.21 


3.81 


-40.13 


5.41 


A 9^ 


A 9^ 


A 11 


M05-2X 


0.56 


-46.76 


3.24 


-31.99 


3.59 


-39.35 


4.63 


3.82 


3.82 


3.87 


cjB97 


0.00—1.00 


A K GO 

-4o. yz 


o /in 


-oo.Do 


o.Zo 


A 1 Rti 
-41.00 


o.y4 


4.86 


4.86 


5.20 


wB97X 


0.16—1.00 


-46.13 


2.61 


-34.34 


5.94 


-42.20 


7.48 


5.34 


5.34 


5.72 


WB97X-D 


0.22—1.00 


-45.71 


2.19 


-35.33 


6.93 


-43.14 


8.42 


5.85 


5.85 


6.42 


wB97X-2(LP) 0.68—1.00 


-45.42 


1.90 


-29.17 


0.77 


-37.83 


3.11 


1.93 


1.93 


2.15 


M06HF 


1.00 


-48.36 


4.84 


-28.23 


-0.17 


-35.75 


1.03 


1.90 


2.01 


2.86 


MP2 


1.00 


-43.95 


0.43 


-25.16 


-3.24 


-30.03 


-4.69 


-2.50 


2.79 


3.30 


CCSD(T) a 


1.00 


-43.52 b 


0.00 


-28.39 


0.00 


-34.72 


0.00 


0.00 


0.00 


0.00 


a The CCSD(T) results, taken from Ref. 

K _ _ . 


0, is 


adopted as the reference. 













' The ZPE corrected binding energy of the proton transferred structure calculated by CCSD(T) in Ref. \1M is inconsistent with 
Ref. However, adopting the geometry of Ref. [l| will yield the results that are consistent with Ref. [12j ]. 



by 



The HF exchange included in the wB97 series is given and 



where 



™B97 scries _ t^LR-HF , n t^SR-HF 
^HF exchange ~ tL "i c i ) 



e HF = 4EE J /cw) 

erf(cjri2) 



'l\2 



(1) 



-V'jffCri)^^)*"!^, (2) 



SR-HF 



erfc(u;ri 2 ) 
ri2 



Cr(riW«r(ra) 



(3) 



Here ri 2 = |ri 2 | = |ri — r 2 | (atomic units are used 
throughout this paper). The parameter ui defines the 
range of the splitting operators. The coefficients for the 
a;B97 series are listed in Table IIIII Since the fraction of 
HF exchange in wB97 series depends on the interelec- 
tronic distance ri2, the trend mentioned previously is 
not as obvious as the global hybrid functionals. But it is 
clear that the wB97X-2(LP), a LC double-hybrid func- 
tional, gives the most accurate results compare to the 
other functionals in the wB97 series. 



As mentioned previously, functionals with large frac- 
tions of HF exchange may perform well for the hemi- 
bonded structure where the serious SIE takes place, they 
may perform unsatisfactorily for the other structures. 
Therefore, we also consider another criterion: the rel- 
ative energies between the three structures, as proposed 
by Cheng et al. QJ] • 



B. Criterion II: Relative energies 

In this criterion, the ground-state energy of the proton 
transferred structure is set to the zero point, i.e. both 
of the ground state energies of the hemibonded structure 
and that of the transition state are relative to the proton 
transferred structure, as shown in Table [TTJ The previ- 
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TABLE II. Relative energies (in kcal/mol) between those three structures of the water dimer radical cation. 



Proton transferred structure Transition state Hemibonded structure 



A jT 4- 1-. ^ J 

Method 


a H F 


hj 


Error 


T 


Error 


T ," 

£j 


Error 


MSE MAE RMS 




0.00 


0.00 


0.00 






-7.27 


16.07 








PBE 


0.00 


0.00 


0.00 






-6.56 


15.36 








M06L 


0.00 


0.00 


0.00 


4.14 


10.86 


-3.24 


12.04 


11 A 
11.40 


11 .40 


11/17 
11.4/ 


B97 


0.19 


0.00 


0.00 


6.93 


8.07 


-0.55 


9.35 


ft 71 
O. 1 1 


ft 71 
o. ( 1 


ft 7 A 
o. ( 4 


B3LYP 


0.20 


0.00 


0.00 


7.51 


7.49 


0.27 


8.53 


s ni 

o.Ul 


ft ni 

o.Ul 


ft 09 

o.Uz 


PBEO 


0.25 


0.00 


0.00 


9.86 


5.14 


2.66 


6.14 


0.04 


^ fiA 
0.04 


0.00 


M06 


0.27 


0.00 


0.00 


9.70 


5.30 


3.40 


5.40 


O.OO 


O.OO 


O.oO 


M05 


0.28 


0.00 


0.00 


9.89 


5.11 


4.04 


4.76 


4.y4 


A OA 

4.y4 


A OA 

4.y4 


TJTTP TTT \/"T> 

BH&HLYP 


0.50 


0.00 


0.00 


16.41 


-1.41 


10.86 


-2.06 


1 7/1 
-1.(4 


1 7 A 
1.(4 


1 77 


B2PLYP 


0.53 


0.00 


0.00 


12.46 


2.54 


4.08 


4.72 


O.OO 


O.OO 


^ 7Q 
o. ( y 


M06-2X 


0.54 


0.00 


0.00 


14.85 


0.15 


6.92 


1.88 


1 09 


1 09 
l.UZ 


1 ^ 

l.OO 


M05-2X 


0.56 


0.00 


0.00 


14.77 


0.23 


7.41 


1.39 


n fti 

U.ol 


O ft1 
U.ol 


1 OO 
l.UU 


cjB97 


0.00—1.00 


0.00 


0.00 


12.29 


2.71 


4.27 


4.53 


3.62 


3.62 


3.73 


cjB97X 


0.16—1.00 


0.00 


0.00 


11.79 


3.21 


3.92 


4.88 


4.05 


4.05 


4.13 


wB97X-D 


0.22—1.00 


0.00 


0.00 


10.38 


4.62 


2.57 


6.23 


5.42 


5.42 


5.48 


cjB97X-2(LP) 0.68—1.00 


0.00 


0.00 


16.25 


-1.25 


7.59 


1.21 


-0.02 


1.23 


1.23 


M06HF 


1.00 


0.00 


0.00 


20.12 


-5.13 


12.61 


-3.81 


-4.47 


4.47 


4.52 


MP 2 


1.00 


0.00 


0.00 


18.79 


-3.79 


13.92 


-5.12 


-4.45 


4.45 


4.50 


CCSD(T) a 


1.00 


0.00 


0.00 


15.14 


0.00 


8.80 


0.00 


0.00 


0.00 


0.00 



a The CCSD(T) results, taken from Ref. [12], is adopted as the reference. 



ously recommended functional, M06HF, performs poorly 
here. 

In this criterion, it is obvious that functionals without 
the exact HF exchange, such as BLYP and PBE, oversta- 
bilize the hemibonded structure and wrongly predict the 
hemibonded structure to be more stable than the proton 
transferred one. 

In addition to the previously recommended function- 
als based on Criterion I, the M05-2X and the M06-2X 
functionals also give accurate relative energies here. Al- 
though they give results that are not accurate enough for 



the binding energies, they yield good relative energies be- 
tween those three structures of the water dimer radical 
cation. 

In fact, functionals that are unable to give reasonable 
binding energies for the hemibonded structure may be 
traced back to the predicted dissociation curves of the 
hemibonded systems. Since many functionals fail to dis- 
sociate the hemibonded structure of water dimer radical 
cation into the correct fragments. The definition of the 
binding energy is not well-defined. Therefore, the entire 
dissociation curve from the hemibonded structure should 
be concerned. 



r 



TABLE III. The coefficients of the SR HF exchange and u> for 
ojB97 series. 



cjB97 



o;B97X cjB97X-D cjB97X-2(LP) 



w(bohr _1 ) 0.4 
C x 0.00 



0.3 
0.16 



0.2 
0.22 



0.3 
0.68 



C. Criterion III: Dissociation behavior 

Due to the severe SIEs associated with DFAs, sys- 
tems with three-electron hemibonds, such as the hemi- 
bonded structure of water dimer radical cation, are espe- 



cially difficult for conventional density functionals. Many 
XC functional cannot dissociate it into the correct frag- 
ments, H2O and H2O" 1 " (ionic state), i.e. they predict that 
the hemibonded structure should be dissociated into two 
fragments, each of which carries half of positive charge 
(covalent state). Fig. [2] shows the dissociation curves 
calculated by various XC functionals. Note that the dis- 
continuous points in Fig. [2] near R = 2.5 angstrom for the 
wB97X-2(LP) and 3 angstrom for the M06HF functional 
are the respective broken-symmetry points. 

The BH&HLYP functional, which has been suggested 
in the equilibrium ground-state energy calculation by the 
earlier reports 0,[Tl| has a spurious barrier on its dissoci- 
ation curve. Although we do not present the dissociation 
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curve of the MPW1K functional, we expect it will suffer 
from the same problem as BH&HLYP. 

The spurious barrier can be removed if the 100% exact 
exchange is adopted in a functional (e.g. M06HF), but its 
shortcoming is described in the previous subsection and is 
thus not recommended. This shortcoming can be greatly 



reduced by the use of LC hybrid functionals, such as 
the wB97 series or the other LC hybrid functionals [33j . 
The LC functionals retain the full HF exchange at the 
long range, while the good cancelation of errors between 
the semilocal exchange and correlation functionals are 
retained at the short range [26[. 




(b) 



20 

10 
0- 
-10- 
-20- 
-30- 



-CCSD(T) 
-MP2 

0)B97X-2(LP) 
- B2PLYP 




123456789 10 

R(angstrom) 

FIG. 2. (a) Dissociation curves for the hemibonded structure 
calculated by various XC functionals. (b) Dissociation curves 
for the hemibonded structure calculated by MP2 and double- 
hybrid functionals. 

In the following, we will explain why LC hybrid func- 
tionals do not suffer from a spurious energy barrier as 
global hybrid functionals do. An estimate of the SIE of 
symmetric radical cations by global hybrid or pure den- 
sity functionals has been derived as: HH 



E 



SIE 



(1 - cvhf) 





i " 








4^. 



(4) 



where C 2" 1 / 3 , (0.5 - C) « -0.29, and J is the 
Coulomb self-interaction energy for the ionic state (in 



the case that the bond electron is localized at either of 
the two fragments). For three-electron-bonded radical 
cations, the bonding between the fragments is accom- 
plished by the delocalized j3 electron, which dominates 
the total SIE [loj]. 

We have derived an estimate of the SIE of symmetric 
radical cations by LC hybrid functionals, with the details 
arranged in the appendix. And the result is 



E 



SIE 



(1 - C x ) 



B(u)C J bH H + 



erfc(wi?) 
~4~R 



( 5 ) 

where B(ui) and J sr (uj) are constants with respect to R. 
But they depend on id, i.e. for LC hybrid functionals with 
different id values, their B(uS) and J SR (u>) are different. 
The dependence of J (id) on ui is defined by 



J (id) 



1 



erfc(wr 12 ) 

., , , PM r i) pp(r2)dndr 2 , (6) 



and for small id, 



B(oj) « 1 - 0.254w(bohr _1 ). 



(7) 



Note that estimate (j4|) is the special case with vanish- 
ing id of estimate (JSJ) - We apply estimate §5§ and ([7} to 
the simplest three-electron-bonded system, and the esti- 
mated HeJ (which is also a three-electron hemibonded 
system) dissociation curves for pure DFT and LC hy- 
brid functional are shown in Fig.[3ja). For simplicity, C x 
is set to zero and the LDA orbital is used for evaluat- 
ing the Coulomb self-interactions; id is set to 0.4 bohr -1 
for the LC hybrid functional. When the exact dissoci- 
ation curve approaches zero, the SIE of the LC hybrid 
functional is already close to a constant, while the SIE 
of pure DFT is still decreasing. Therefore the dissoci- 
ation curve by the LC hybrid functional does not dis- 
play a spurious energy barrier as that of the pure DFT. 
Another effect of the long-range correction is the reduc- 
tion of Coulomb self-interaction for the ionic state. In 
Fig. EJa), (0.5 — C)J w —92 kcal/mol has been modified 
to [0.5 - B(id)C] J SR (u>) « -40 kcal/mol. 

The formal feature of spurious energy barrier is one 
more turning point (at which the derivative changes sign) 
on top of the barrier, in addition to the one in the po- 
tential well. Since the ground state energy calculated 
by a functional is approximately the exact ground state 
energy plus the SIE produced by that functional, 



E 



DFT 



E c 



E 



SIE 



(8) 
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turning points occur when |<iE SIE /(iR| equals the deriva- 
tive of the exact curve, i.e. points where \dE sm / dR\ in- 
tersects the derivative of the exact curve. \dE slE /dR\ by 
pure DFT is l/4i? 2 . For the LC hybrid functional, there 
is an extra multiplicative factor: 



(a) 



4iT 



dE 



SIE 



dR 



2ujR 



exp [-(wi?) 2 ] + erfc(wi?). (9) 



As shown in Fig. [3jb) , for sufficient u> value (0.2 bohr -1 
or more), this factor can change the decay nature of the 
derivative magnitude, from power law to exponential. 
Thus, the SIE derivative magnitude curve of typical LC 
hybrid functionals can avoid second intersection with the 
derivative of the exact curve. Global hybrid functionals 
simply scale down the SIE derivative curve by a constant, 
so they cannot avoid the second intersection, unless ohf 
approaches unity. This explains why global hybrid func- 
tionals display a spurious energy barrier which LC hybrid 
functionals avoid. 

Note that LC hybrid functionals which do not contain 
SR HF exchange, such as wB97, still can be free from 
the spurious barrier, but will lose the possibility to pre- 
dict symmetry breaking during the dissociation, i.e. the 
dissociation curve can not converge to zero. This means 
that the SR HF exchange is also important. In fact, 
the covalent (symmetric) state and the ionic (symmetry- 
broken) state are nearly degenerate by CCSD(T) calcula- 
tions [lfj . But most of the XC functionals cannot predict 
that these two states are degenerate: Due to the serious 
SIE, they usually overstabilize the covalent state. There- 
fore, functional which can predict that the ionic state 
is more stable than the covalent state (i.e. the hemi- 
bonded structure will dissociate into H 2 and H 2 + ) will 
give correct dissociation limit. Very recently, a double- 
hybrid functional containing a very large fraction of HF 
exchange (sa 79%) [36|, has been shown to be promising 
for reducing the SIEs in hemibonded systems. 

The dissociation curves of the hemibonded structure 
calculated by double-hybrid functionals and MP2 are 
shown in Fig.[2Jb). Note that the PT2 calculation should 
be executed in the stable wave function. For example, al- 
though the dissociation curve of the covalent state seems 
more stable than that of the ionic state, we should choose 
the dissociation curve of the ionic state as an actual dis- 
sociation behavior calculated by MP2. Since the HF the- 
ory, which provides the reference orbitals for computing 
the MP 2 correlation energy, predicts the ionic state to 
be more stable than the covalent one. Thus we choose 
the dissociation curve of the ionic state rather than that 
of the covalent state. Fig. |2Jb) shows that the wB97X- 
2 (LP) functional and the MP2 theory can predict the 
correct dissociation limits. While the B2PLYP functional 
cannot. 

There is another way to define the dissociation behav- 
ior of the XC functionals: Since many of the XC function- 
als cannot predict that the hemibonded structure of the 
ionized water dimer will dissociate into H2O and H2 + , 
we set the zeros of dissociation curves to their respective 




2345 
R (angstrom) 



(b) 




CCSD(T) 

- - |dSIE(LC)/dR| 
|dSIE(pure)/dR| 



R (angstrom) 



FIG. 3. The spurious energy barrier of the hemibonded sys- 
tems predicted by DFT functionals can be illustrated quali- 
tatively, (a) Comparison of HeJ dissociation curves by pure 
DFT and LC hybrid functional using estimate ([5}. Zero level 
is set to E(He)+E(HeJ) for each method, (b) First derivative 
of exact HeJ dissociation curve and SIE derivative magni- 
tude curves by pure DFT and LC hybrid functional using 
estimate 



dissociation limits, as shown in Fig. 01 In this definition, 
we focus on the potential curve experienced by the two 
fragments during the dissociation process of the hemi- 
bonded structure. Surprisingly, the dissociation curve 
of the wB97 functional is extremely close to that of the 
CCSD(T) theory. This means that the wB97 gives the 
best potential energy curve toward the dissociation pro- 
cess. The previous suggested functionals [si, [Hj], such as 
BH&HLYP, yield potential curves that are too shallow. 
Functionals which predict symmetry-breaking solutions 
during the dissociation process (e.g. M06HF and wB97X- 
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2(LP)) are found to yield the dissociation curves that are 
narrower than that of the CCSD(T) theory. 




4 5 6 

R(angstrom) 



FIG. 4. Dissociation curves for the hemibonded structure. 
The zeros of dissociation curves are set to their respective 
dissociation limits. 

Finally, we discuss the dissociation of the proton trans- 
ferred structure of the water dimer cation. In this struc- 
ture, the SIEs associated with functionals for this struc- 
ture are not as large as those for the hemibonded struc- 
ture, so all the dissociation curves are very similar, as 
shown in Fig. [SJ 




R(angstrom) 



FIG. 5. The dissociation curves for the proton transferred 
structure of ionized water dimer. 

The results of the water dimer radical cation using 
three criteria we proposed are summarized in Table. IIVI 
We find that the functional which performs well based on 
these three criteria is the u;B97X-2(LP) functional, yield- 
ing the accurate binding energies, relative energies, and 
the correct dissociation limit. However, this functional 
yields a dissociation curve of the hemibonded structure 
that is a little narrower than that of CCSD(T). For ap- 
plications sensitive to the shape of potential of the hemi- 
bonded structure, we suggest to use wB97: Although this 



functional neither gives a dissociation curve of the hemi- 
bonded structure that converges to zero nor yields accu- 
rate binding energies, this functional gives a dissociation 
curve which has nearly the same shape as the curve calcu- 
lated by CCSD(T). Thus, we recommend this functional 
for researchers who like to perform the molecular dynam- 
ics of the water dimer cation. 



IV. CONCLUSIONS 

We have proposed three criteria to examine the per- 
formance of density functionals on water dimer radical 
cation, and explained why LC hybrid functionals gener- 
ally work better than conventional density functionals for 
hemibonded systems. 

The previously recommended functional, BH&HLYP, 
cannot dissociate the hemibonded structure of the water 
dimer cation into the correct fragments: H2O and H2 + . 
Furthermore, the BH&HLYP dissociation curve displays 
an unphysical repulsive barrier, and is too shallow for 
molecular dynamics simulations. Such a spurious barrier 
could be removed by functionals with very large fractions 
of HF exchange (e.g. M06HF and wB97X-2(LP)). LC hy- 
brid functionals, such as the wB97 series, are shown to be 
accurate for dissociation curves of the hemibonded struc- 
ture (i.e. no spurious barriers), and thus are suitable for 
molecular dynamics simulations of larger-size systems. 
For research which is sensitive to the dissociation curve 
experienced by the fragments of the hemibonded struc- 
ture, we recommend the use of the wB97 functional. For 
researchers who do not sure which criterion is the most 
important factor during the simulation of water dimer 
radical cation, we recommend the use of the wB97X- 
2 (LP) functional, which is overall good throughout those 
three criteria we proposed. 



APPENDIX 

We first reproduce estimate ((H) for the SIE of sym- 
metric radical cations by global hybrid or pure density 
functionals HH . Since most XC functionals predict 
that the covalent state is more stable than the ionic state, 
the SIE of the covalent state is of interest. Neglecting the 
small SIE in the correlation energy, we have 



pSIE ^ to 
cov ~ ^cov 



E: 



SI 



Because the SIE is small in the ionic state, 



E 



si 



x. ionic 



J 



SI 



(Al) 



(A2) 



it is favorable to express J^ v and E^ 1 cov in terms of 
the Coulomb self-interaction energy for the ionic state 
Jfonic> or simply J in subsection IIII CI This can be done 
by substituting the density of the delocalized (3 electron 
(which causes the serious SIE of the hemibonding systems 
[lfj) in the covalent state with the density of that electron 
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in the ionic state. Note that there are two situations for 
the ionic state: one is the electric charge localized in 
fragment A and the other is the electric charge localized 
in fragment B. To a good approximation the density of 
the delocalized (3 electron of the covalent state can be 
expressed as 



rcov 



(r) 



(r) 



(A3) 



The Coulomb self-interaction for the covalent state can 
be expressed as 



1 



rSI 



1 



P 1 sa - TP 

cov 2 ionic 1 ^^j> i 



(A4) 



if one assumes that R is large compared to the spatial 
extent of and pg. This expression can be applied to 
the self-interaction HF exchange energy for the covalent 
state, 



E; 



SI, HF 



J. 



SI 
cov ' 



(A5) 



The pure-DFT self-exchange energy for the covalent 
state is 



For LDA, C = 2" 1 / 3 w 0.79. Combining 

Ac, cov - tt HFfi x . cov + I 1 _ "HFj-C'x, cov 



(A6) 



(A7) 



and estimate (|A4j) . one can obtain estimate ((4]). 

The SIE estimate for LC hybrid functionals can be 
derived from the same manner as the above one for global 
hybrid functionals. Substituting estimate (|A3[) into the 
self-interaction LR HF exchange yields 



R 



SI. LR— HF 



:•/; 



SI, LR 



erf(wR) 
AR 



(A8) 



where we define 

7 si, SR/ \ 1 ff i ,erfc(o;ri2) 
^ionic ( w ) = 2 J J P ^ Y ^> ~ 2 M r 2)dridr 2 . 

(All) 

The SR-DFA self-exchange energy for the covalent 
state is 



pSI, SR-DFA 
v x, cov . 

For SR-LDA 126, 3 



si. 



(A12) 



(A13) 

kp/3 = (67r 2 p j a(r)) 1 / 3 is the local Fermi wave vector, and 
the attenuation function is given by 



9 \ 

F(A) = l-^[2V^erf(A- 1 ) + 

i 4 V5r t n \ 

k1 , for small A. 



(A14) 



For small u, and with the density approximated as one 
electron in a sphere with Bohr radius ao, 

B(u) « 1 - 3 _B / 3 4^2(^2- l)^Fwao 

w 1 - 0.254w(bohr" 1 ). (A15) 

Combining 



E =E 

x. cov X. cov 



SI. LR-HF 



SI, SR-HF 



+ (1 - C x )El 



SI, SR-DFA 
cov 



(A16) 



and estimate (|A4j) . one can obtain estimate ([5]). 



where we define 



tSL LR/ n 1 
J ionic M = 2 



^( r i) erf(U,? ' 12) ^( r 2)rfr 1 rfr 2 . (A9) 



ri2 



Since the integration of ^nic^ '( w ) is on ly over one frag- 
ment, it is independent of R. Likewise, the SR HF ex- 
change of the covalent state is 



E. 



SI, SR-HF 



1 r Si, SR/ a , erfc(wi?) 

2 J ionic l W J H ^ 



(A10) 
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TABLE IV. Summary of results based on the three criteria. 
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a The CCSD(T) results, taken from Ref. [12|. is adopted as the reference. 



